Near Analysis Models



What Was Tested: Do crop yields increase from GCOs globally?



What This Allows Me To Do: Identify if yield-distance relationships (near analysis) are a function of evolutionary escape or other variables.



What Are The Model Assumptions: GLS models are used as they have less assumptions about error structure and heteroscadicity.

#Barley -

#Select Data and Rescale Distance
near<-read.csv(file = "Near_All/Barley_Near_Centroid.csv", header=T)
near$rescale_ND <- near$NEAR_DIST/(1000*1000)
barley.near<-near %>% select(mean_barle,rescale_ND)
barley.near<-barley.near %>% filter(mean_barle > 0, na.rm=TRUE)
barley.near<-barley.near %>% filter(rescale_ND > 0, na.rm=TRUE)
barley.near$logBarley<-log10(barley.near$mean_barle)

#Plot
#Plot
ggplot(barley.near, aes(x=rescale_ND, y=logBarley)) +
  geom_point(alpha=0.1, color="#00A9B4") + geom_smooth(method = "lm",color="#DB3500") + labs(color='GCO') +theme_justin +xlab("Distance From GCO (1000's km)") +ylab ("Log10 Yield")
## `geom_smooth()` using formula 'y ~ x'
Barley Near near Model,

Barley Near near Model, 

###  Model - Barley 
# regular OLS no variance structure
barley.ols <- gls(logBarley ~ rescale_ND, data = barley.near)
# varFixed (variance changes linearly with X)
barley.fixed <- update(barley.ols, .~., weights = varFixed(~rescale_ND))
# varPower (variance changes as a power function with X)
barley.power <- update(barley.ols, . ~ ., weights = varPower(form = ~rescale_ND))
# varExp (variance changes as an exponential function of x)
barley.exp <- update(barley.ols, . ~., weights = varExp(form = ~ rescale_ND)) # error
# varConstPower (constant plus a power function of X (useful if X includes 0))
barley.ConstPower <- update(barley.ols, . ~., weights = varConstPower(form = ~ rescale_ND))
# compare all models by AIC
AIC(barley.ols, barley.fixed, barley.power, barley.ConstPower, barley.exp)
##                   df      AIC
## barley.ols         3 19822.15
## barley.fixed       3 25277.66
## barley.power       4 19729.77
## barley.ConstPower  5 19731.77
## barley.exp         4 19730.85
#varPower
summary(barley.exp)
## Generalized least squares fit by REML
##   Model: logBarley ~ rescale_ND 
##   Data: barley.near 
##        AIC      BIC    logLik
##   19730.85 19760.14 -9861.423
## 
## Variance function:
##  Structure: Exponential of variance covariate
##  Formula: ~rescale_ND 
##  Parameter estimates:
##       expon 
## -0.01694584 
## 
## Coefficients:
##                  Value  Std.Error   t-value p-value
## (Intercept) -0.1452756 0.01100551 -13.20026       0
## rescale_ND   0.0147758 0.00138106  10.69889       0
## 
##  Correlation: 
##            (Intr)
## rescale_ND -0.867
## 
## Standardized residuals:
##        Min         Q1        Med         Q3        Max 
## -4.5651801 -0.4496173  0.2959025  0.7252327  1.5929104 
## 
## Residual standard error: 0.6499998 
## Degrees of freedom: 11197 total; 11195 residual

#Cassava -

#Select Data and Rescale Distance
near<-read.csv(file = "Near_All/Cassava_Near_Centroid.csv", header=T)
near$rescale_ND <- near$NEAR_DIST/(1000*1000)
cassava.near<-near %>% select(mean_cassa,rescale_ND)
cassava.near<-cassava.near %>% filter(mean_cassa > 0, na.rm=TRUE)
cassava.near<-cassava.near %>% filter(rescale_ND > 0, na.rm=TRUE)
cassava.near$logCassava<-log10(cassava.near$mean_cassa)

summary(cassava.near)
##    mean_cassa         rescale_ND         logCassava     
##  Min.   : 0.01222   Min.   : 0.02601   Min.   :-1.9128  
##  1st Qu.: 3.67075   1st Qu.: 3.44972   1st Qu.: 0.5648  
##  Median : 7.73612   Median :10.39709   Median : 0.8885  
##  Mean   : 8.17809   Mean   :10.58149   Mean   : 0.7396  
##  3rd Qu.:12.24387   3rd Qu.:16.97533   3rd Qu.: 1.0879  
##  Max.   :38.89642   Max.   :19.72692   Max.   : 1.5899
#Plot
ggplot(cassava.near, aes(x=rescale_ND, y=logCassava)) +
  geom_point(alpha=0.2, color="#00A9B4") + geom_smooth(method = "lm",color="#DB3500") +theme_justin +xlab("Distance From GCO (1000's km)") +ylab ("Log10 Yield")
## `geom_smooth()` using formula 'y ~ x'
Cassava Near near Model,

Cassava Near near Model, 

###  Model - Cassava 
# regular OLS no variance structure
cassava.ols <- gls(logCassava ~ rescale_ND, data = cassava.near)
# varFixed (variance changes linearly with X)
cassava.fixed <- update(cassava.ols, .~., weights = varFixed(~rescale_ND))
# varPower (variance changes as a power function with X)
cassava.power <- update(cassava.ols, . ~ ., weights = varPower(form = ~rescale_ND))
# varExp (variance changes as an exponential function of x)
cassava.exp <- update(cassava.ols, . ~., weights = varExp(form = ~ rescale_ND)) # error
# varConstPower (constant plus a power function of X (useful if X includes 0))
cassava.ConstPower <- update(cassava.ols, . ~., weights = varConstPower(form = ~ rescale_ND))
# compare all models by AIC
AIC(cassava.ols, cassava.fixed, cassava.power, cassava.ConstPower, cassava.exp)
##                    df       AIC
## cassava.ols         3  8900.148
## cassava.fixed       3 13094.247
## cassava.power       4  8862.351
## cassava.ConstPower  5  8864.351
## cassava.exp         4  8756.726
#varExp
summary(cassava.exp)
## Generalized least squares fit by REML
##   Model: logCassava ~ rescale_ND 
##   Data: cassava.near 
##        AIC      BIC    logLik
##   8756.726 8783.727 -4374.363
## 
## Variance function:
##  Structure: Exponential of variance covariate
##  Formula: ~rescale_ND 
##  Parameter estimates:
##       expon 
## -0.02047272 
## 
## Coefficients:
##                 Value   Std.Error  t-value p-value
## (Intercept) 0.5155451 0.013402171 38.46728       0
## rescale_ND  0.0207906 0.000983958 21.12954       0
## 
##  Correlation: 
##            (Intr)
## rescale_ND -0.895
## 
## Standardized residuals:
##        Min         Q1        Med         Q3        Max 
## -5.5684850 -0.3568539  0.3005624  0.7105753  1.7183862 
## 
## Residual standard error: 0.5998532 
## Degrees of freedom: 6315 total; 6313 residual

#Groundnut -

#Select Data and Rescale Distance
near<-read.csv(file = "Near_All/Groundnut_Near_Centroid.csv", header=T)
near$rescale_ND <- near$NEAR_DIST/(1000*1000)
groundnut.near<-near %>% select(mean_groun,rescale_ND)
groundnut.near<-groundnut.near %>% filter(mean_groun > 0, na.rm=TRUE)
groundnut.near<-groundnut.near %>% filter(rescale_ND > 0, na.rm=TRUE)
groundnut.near$logGroundnut<-log10(groundnut.near$mean_groun)

#Plot
ggplot(groundnut.near, aes(x=rescale_ND, y=logGroundnut)) +
  geom_point(alpha=0.2, color="#00A9B4") + geom_smooth(method = "lm",color="#DB3500") +theme_justin +xlab("Distance From GCO (1000's km)") +ylab ("Log10 Yield")
## `geom_smooth()` using formula 'y ~ x'
Groundnut Near near Model,

Groundnut Near near Model, 

###  Model - Groundnut 
# regular OLS no variance structure
groundnut.ols <- gls(logGroundnut ~ rescale_ND, data = groundnut.near)
# varFixed (variance changes linearly with X)
groundnut.fixed <- update(groundnut.ols, .~., weights = varFixed(~rescale_ND))
# varPower (variance changes as a power function with X)
groundnut.power <- update(groundnut.ols, . ~ ., weights = varPower(form = ~rescale_ND))
# varExp (variance changes as an exponential function of x)
groundnut.exp <- update(groundnut.ols, . ~., weights = varExp(form = ~ rescale_ND)) # error
# varConstPower (constant plus a power function of X (useful if X includes 0))
groundnut.ConstPower <- update(groundnut.ols, . ~., weights = varConstPower(form = ~ rescale_ND))
# compare all models by AIC
AIC(groundnut.ols, groundnut.ConstPower,groundnut.power, groundnut.fixed, groundnut.exp)
##                      df      AIC
## groundnut.ols         3 13520.01
## groundnut.ConstPower  5 13523.68
## groundnut.power       4 13521.68
## groundnut.fixed       3 18936.25
## groundnut.exp         4 13520.79
#varPower
summary(groundnut.exp)
## Generalized least squares fit by REML
##   Model: logGroundnut ~ rescale_ND 
##   Data: groundnut.near 
##        AIC     BIC    logLik
##   13520.79 13548.9 -6756.397
## 
## Variance function:
##  Structure: Exponential of variance covariate
##  Formula: ~rescale_ND 
##  Parameter estimates:
##       expon 
## 0.001516301 
## 
## Coefficients:
##                   Value   Std.Error   t-value p-value
## (Intercept) -0.24712201 0.013296847 -18.58501   0e+00
## rescale_ND   0.00429596 0.001069283   4.01761   1e-04
## 
##  Correlation: 
##            (Intr)
## rescale_ND -0.894
## 
## Standardized residuals:
##        Min         Q1        Med         Q3        Max 
## -4.4934597 -0.4081239  0.2516889  0.6904358  1.8240400 
## 
## Residual standard error: 0.5352301 
## Degrees of freedom: 8321 total; 8319 residual

#Maize -

#Select Data and Rescale Distance
near<-read.csv(file = "Near_All/Maize_Near_Centroid.csv", header=T)
near$rescale_ND <- near$NEAR_DIST/(1000*1000)
maize.near<-near %>% select(mean_maize,rescale_ND)
maize.near<-maize.near %>% filter(mean_maize > 0, na.rm=TRUE)
maize.near<-maize.near %>% filter(rescale_ND > 0, na.rm=TRUE)
maize.near$logMaize<-log10(maize.near$mean_maize)

#Plot
ggplot(maize.near, aes(x=rescale_ND, y=logMaize)) +geom_point(alpha=0.2, color="#F0AB00") +geom_smooth(method = "lm",color="#DB3500") +theme_justin +xlab("Distance From GCO (1000's km)") +ylab ("Log10 Yield")
## `geom_smooth()` using formula 'y ~ x'
Maize Near near Model,

Maize Near near Model, 

###  Model - Maize 
# regular OLS no variance structure
maize.ols <- gls(logMaize ~ rescale_ND, data = maize.near)
# varFixed (variance changes linearly with X)
maize.fixed <- update(maize.ols, .~., weights = varFixed(~rescale_ND))
# varPower (variance changes as a power function with X)
maize.power <- update(maize.ols, . ~ ., weights = varPower(form = ~rescale_ND))
# varExp (variance changes as an exponential function of x) - Does not Converge
maize.exp <- update(maize.ols, . ~., weights = varExp(form = ~ rescale_ND)) # error
# varConstPower (constant plus a power function of X (useful if X includes 0))
maize.ConstPower <- update(maize.ols, . ~., weights = varConstPower(form = ~ rescale_ND))
# compare all models by AIC
AIC(maize.ols, maize.fixed, maize.power, maize.ConstPower,maize.exp)
##                  df      AIC
## maize.ols         3 23829.13
## maize.fixed       3 26600.48
## maize.power       4 23786.72
## maize.ConstPower  5 23788.72
## maize.exp         4 23830.41
#varExp
summary(maize.exp)
## Generalized least squares fit by REML
##   Model: logMaize ~ rescale_ND 
##   Data: maize.near 
##        AIC      BIC   logLik
##   23830.41 23860.54 -11911.2
## 
## Variance function:
##  Structure: Exponential of variance covariate
##  Formula: ~rescale_ND 
##  Parameter estimates:
##       expon 
## 0.001294945 
## 
## Coefficients:
##                  Value  Std.Error   t-value p-value
## (Intercept)  0.4023974 0.01228996  32.74195       0
## rescale_ND  -0.0264575 0.00112441 -23.53007       0
## 
##  Correlation: 
##            (Intr)
## rescale_ND -0.918
## 
## Standardized residuals:
##        Min         Q1        Med         Q3        Max 
## -4.6574695 -0.4741918  0.2375453  0.7011139  2.0342761 
## 
## Residual standard error: 0.5660735 
## Degrees of freedom: 13792 total; 13790 residual

#Millet -

#Select Data and Rescale Distance
near<-read.csv(file="Near_All/Millet_Near_Centroid.csv", header=T)
near$rescale_ND <- near$NEAR_DIST/(1000*1000)
millet.near<-near %>% select(mean_mille,rescale_ND)
millet.near<-millet.near %>% filter(mean_mille > 0, na.rm=TRUE)
millet.near<-millet.near %>% filter(rescale_ND > 0, na.rm=TRUE)
millet.near$logMillet<-log10(millet.near$mean_mille)

#Plot
ggplot(millet.near, aes(x=rescale_ND, y=logMillet)) + geom_point(alpha=0.2, color="#FFB100") + geom_smooth(method = "lm",color="#DB3500") +theme_justin +xlab("Distance From GCO (1000's km)") +ylab ("Log10 Yield")
## `geom_smooth()` using formula 'y ~ x'
Millet Near near Model,

Millet Near near Model, 

###  Model - Millet 
# regular OLS no variance structure
millet.ols <- gls(logMillet ~ rescale_ND, data = millet.near)
# varFixed (variance changes linearly with X)
millet.fixed <- update(millet.ols, .~., weights = varFixed(~rescale_ND))
# varPower (variance changes as a power function with X)
millet.power <- update(millet.ols, . ~ ., weights = varPower(form = ~rescale_ND))
# varExp (variance changes as an exponential function of x)
millet.exp <- update(millet.ols, . ~., weights = varExp(form = ~ rescale_ND)) # error
# varConstPower (constant plus a power function of X (useful if X includes 0))
millet.ConstPower <- update(millet.ols, . ~., weights = varConstPower(form = ~ rescale_ND))
# compare all models by AIC
AIC(millet.ols, millet.fixed, millet.power, millet.ConstPower, millet.exp)
##                   df      AIC
## millet.ols         3 12847.29
## millet.fixed       3 15549.03
## millet.power       4 12843.80
## millet.ConstPower  5 12845.80
## millet.exp         4 12849.22
#varPower
summary(millet.power)
## Generalized least squares fit by REML
##   Model: logMillet ~ rescale_ND 
##   Data: millet.near 
##       AIC      BIC    logLik
##   12843.8 12871.38 -6417.901
## 
## Variance function:
##  Structure: Power of variance covariate
##  Formula: ~rescale_ND 
##  Parameter estimates:
##      power 
## 0.02436811 
## 
## Coefficients:
##                  Value   Std.Error    t-value p-value
## (Intercept) -0.3207978 0.011790772 -27.207533   0e+00
## rescale_ND  -0.0060439 0.001533657  -3.940818   1e-04
## 
##  Correlation: 
##            (Intr)
## rescale_ND -0.816
## 
## Standardized residuals:
##        Min         Q1        Med         Q3        Max 
## -4.9859199 -0.4021284  0.2580831  0.6758465  1.7350781 
## 
## Residual standard error: 0.560494 
## Degrees of freedom: 7298 total; 7296 residual

#Rapeseed -

#Select Data and Rescale Distance
near<-read.csv(file="./Near_All/Rapeseed_Near_Centroid.csv", header=T)
near$rescale_ND <- near$NEAR_DIST/(1000*1000)
rapeseed.near<-near %>% select(mean_rapes,rescale_ND)
rapeseed.near<-rapeseed.near %>% filter(mean_rapes > 0, na.rm=TRUE)
rapeseed.near<-rapeseed.near %>% filter(rescale_ND > 0, na.rm=TRUE)
rapeseed.near$logRapeseed<-log10(rapeseed.near$mean_rapes)

#Plot
ggplot(rapeseed.near, aes(x=rescale_ND, y=logRapeseed)) +geom_point(alpha=0.2, color="#FFAF02") + geom_smooth(method = "lm",color="#DB3500") +theme_justin +xlab("Distance From GCO (1000's km)") +ylab ("Log10 Yield")
## `geom_smooth()` using formula 'y ~ x'
Rapeseed Near near Model,

Rapeseed Near near Model, 

###  Model - Rapeseed 
# regular OLS no variance structure
rapeseed.ols <- gls(logRapeseed ~ rescale_ND, data = rapeseed.near)
# varFixed (variance changes linearly with X)
rapeseed.fixed <- update(rapeseed.ols, .~., weights = varFixed(~rescale_ND))
# varPower (variance changes as a power function with X)
rapeseed.power <- update(rapeseed.ols, . ~ ., weights = varPower(form = ~rescale_ND))
# varExp (variance changes as an exponential function of x)
rapeseed.exp <- update(rapeseed.ols, . ~., weights = varExp(form = ~ rescale_ND)) # error
# varConstPower (constant plus a power function of X (useful if X includes 0))
rapeseed.ConstPower <- update(rapeseed.ols, . ~., weights = varConstPower(form = ~ rescale_ND))
# compare all models by AIC
AIC(rapeseed.ols, rapeseed.fixed, rapeseed.power, rapeseed.ConstPower, rapeseed.exp)
##                     df      AIC
## rapeseed.ols         3 15008.97
## rapeseed.fixed       3 19957.49
## rapeseed.power       4 15008.34
## rapeseed.ConstPower  5 15010.34
## rapeseed.exp         4 15003.10
#varPower
summary(rapeseed.exp)
## Generalized least squares fit by REML
##   Model: logRapeseed ~ rescale_ND 
##   Data: rapeseed.near 
##       AIC      BIC    logLik
##   15003.1 15031.64 -7497.552
## 
## Variance function:
##  Structure: Exponential of variance covariate
##  Formula: ~rescale_ND 
##  Parameter estimates:
##        expon 
## -0.005850967 
## 
## Coefficients:
##                   Value   Std.Error    t-value p-value
## (Intercept) -0.24958856 0.010965742 -22.760754   0e+00
## rescale_ND  -0.00563787 0.001513497  -3.725062   2e-04
## 
##  Correlation: 
##            (Intr)
## rescale_ND -0.857
## 
## Standardized residuals:
##        Min         Q1        Med         Q3        Max 
## -4.5076091 -0.4305570  0.3144485  0.7080704  1.5508114 
## 
## Residual standard error: 0.5629072 
## Degrees of freedom: 9258 total; 9256 residual

#Rice -

#Select Data and Rescale Distance
near<-read.csv(file="./Near_All/Rice_Near_Centroid.csv", header=T)
near$rescale_ND <- near$NEAR_DIST/(1000*1000)
rice.near<-near %>% select(mean_rice,rescale_ND)
rice.near<-rice.near %>% filter(mean_rice > 0, na.rm=TRUE)
rice.near<-rice.near %>% filter(rescale_ND > 0, na.rm=TRUE)
rice.near$logRice<-log10(rice.near$mean_rice)

#Plot
ggplot(rice.near, aes(x=rescale_ND, y=logRice)) +geom_point(alpha=0.2, color="#FBA900") + geom_smooth(method = "lm",color="#DB3500") +theme_justin +xlab("Distance From GCO (1000's km)") +ylab ("Log10 Yield")
## `geom_smooth()` using formula 'y ~ x'
Rice Near near Model,

Rice Near near Model, 

###  Model - Rice 
# regular OLS no variance structure
rice.ols <- gls(logRice ~ rescale_ND, data = rice.near)
# varFixed (variance changes linearly with X)
rice.fixed <- update(rice.ols, .~., weights = varFixed(~rescale_ND))
# varPower (variance changes as a power function with X)
rice.power <- update(rice.ols, . ~ ., weights = varPower(form = ~rescale_ND))
# varExp (variance changes as an exponential function of x)
rice.exp <- update(rice.ols, . ~., weights = varExp(form = ~ rescale_ND)) # error
# varConstPower (constant plus a power function of X (useful if X includes 0))
rice.ConstPower <- update(rice.ols, . ~., weights = varConstPower(form = ~ rescale_ND))
# compare all models by AIC
AIC(rice.ols, rice.fixed, rice.power, rice.ConstPower, rice.exp)
##                 df      AIC
## rice.ols         3 13556.55
## rice.fixed       3 14019.13
## rice.power       4 13242.15
## rice.ConstPower  5 13244.15
## rice.exp         4 13449.33
#varPower
summary(rice.power)
## Generalized least squares fit by REML
##   Model: logRice ~ rescale_ND 
##   Data: rice.near 
##        AIC      BIC    logLik
##   13242.15 13270.31 -6617.077
## 
## Variance function:
##  Structure: Power of variance covariate
##  Formula: ~rescale_ND 
##  Parameter estimates:
##     power 
## 0.1988213 
## 
## Coefficients:
##                  Value   Std.Error   t-value p-value
## (Intercept)  0.3714575 0.008714763  42.62394       0
## rescale_ND  -0.0291549 0.001030728 -28.28577       0
## 
##  Correlation: 
##            (Intr)
## rescale_ND -0.768
## 
## Standardized residuals:
##        Min         Q1        Med         Q3        Max 
## -4.5128640 -0.4174723  0.1789334  0.6418429  1.7784249 
## 
## Residual standard error: 0.3694383 
## Degrees of freedom: 8428 total; 8426 residual

#Rye -

#Select Data and Rescale Distance
near<-read.csv(file="./Near_All/Rye_Near_Centroid.csv", header=T)
near$rescale_ND <- near$NEAR_DIST/(1000*1000)
rye.near<-near %>% select(mean_rye_Y,rescale_ND)
rye.near<-rye.near %>% filter(mean_rye_Y > 0, na.rm=TRUE)
rye.near<-rye.near %>% filter(rescale_ND > 0, na.rm=TRUE)
rye.near$logRye<-log10(rye.near$mean_rye_Y)

#Plot
ggplot(rye.near, aes(x=rescale_ND, y=logRye)) +geom_point(alpha=0.2, color="#FFB200") +geom_smooth(method = "lm",color="#DB3500")  +theme_justin +xlab("Distance From GCO (1000's km)") +ylab ("Log10 Yield")
## `geom_smooth()` using formula 'y ~ x'
Rye Near near Model,

Rye Near near Model, 

###  Model - Rye 
# regular OLS no variance structure
rye.ols <- gls(logRye ~ rescale_ND, data = rye.near)
# varFixed (variance changes linearly with X)
rye.fixed <- update(rye.ols, .~., weights = varFixed(~rescale_ND))
# varPower (variance changes as a power function with X)
rye.power <- update(rye.ols, . ~ ., weights = varPower(form = ~rescale_ND))
# varExp (variance changes as an exponential function of x)
rye.exp <- update(rye.ols, . ~., weights = varExp(form = ~ rescale_ND)) # error
# varConstPower (constant plus a power function of X (useful if X includes 0))
rye.ConstPower <- update(rye.ols, . ~., weights = varConstPower(form = ~ rescale_ND))
# compare all models by AIC
AIC(rye.ols, rye.fixed, rye.power, rye.ConstPower, rye.exp)
##                df      AIC
## rye.ols         3 14495.61
## rye.fixed       3 18866.35
## rye.power       4 14483.81
## rye.ConstPower  5 14367.69
## rye.exp         4 14497.54
#varConstPower
summary(rye.ConstPower)
## Generalized least squares fit by REML
##   Model: logRye ~ rescale_ND 
##   Data: rye.near 
##        AIC      BIC    logLik
##   14367.69 14402.87 -7178.846
## 
## Variance function:
##  Structure: Constant plus power of variance covariate
##  Formula: ~rescale_ND 
##  Parameter estimates:
##     const     power 
##  5.970222 -2.338874 
## 
## Coefficients:
##                   Value   Std.Error    t-value p-value
## (Intercept)  0.02234046 0.010606468   2.106306  0.0352
## rescale_ND  -0.03242745 0.001418423 -22.861619  0.0000
## 
##  Correlation: 
##            (Intr)
## rescale_ND -0.814
## 
## Standardized residuals:
##        Min         Q1        Med         Q3        Max 
## -4.8648858 -0.5176796  0.2805052  0.6733118  2.1347298 
## 
## Residual standard error: 0.0924272 
## Degrees of freedom: 8400 total; 8398 residual

#Sorghum -

#Select Data and Rescale Distance
near<-read.csv(file="./Near_All/Sorghum_Near_Centroid.csv", header=T)
near$rescale_ND <- near$NEAR_DIST/(1000*1000)
sorghum.near<-near %>% select(mean_sorgh,rescale_ND)
sorghum.near<-sorghum.near %>% filter(rescale_ND > 0, na.rm=TRUE)
sorghum.near<-sorghum.near %>% filter(mean_sorgh > 0, na.rm=TRUE)
sorghum.near$logSorghum<-log10(sorghum.near$mean_sorgh)

#Plot
ggplot(sorghum.near, aes(x=rescale_ND, y=logSorghum)) + geom_point(alpha=0.2, color="#00A9B4") +geom_smooth(method = "lm",color="#DB3500") +theme_justin +xlab("Distance From GCO (1000's km)") +ylab ("Log10 Yield")
## `geom_smooth()` using formula 'y ~ x'
Sorghum Near near Model,

Sorghum Near near Model, 

###  Model - Sorghum 
# regular OLS no variance structure
sorghum.ols <- gls(logSorghum ~ rescale_ND, data = sorghum.near)
# varFixed (variance changes linearly with X)
sorghum.fixed <- update(sorghum.ols, .~., weights = varFixed(~rescale_ND))
# varPower (variance changes as a power function with X)
sorghum.power <- update(sorghum.ols, . ~ ., weights = varPower(form = ~rescale_ND))
# varExp (variance changes as an exponential function of x)
sorghum.exp <- update(sorghum.ols, . ~., weights = varExp(form = ~ rescale_ND)) # error
# varConstPower (constant plus a power function of X (useful if X includes 0))
sorghum.ConstPower <- update(sorghum.ols, . ~., weights = varConstPower(form = ~ rescale_ND))
# compare all models by AIC
AIC(sorghum.ols, sorghum.fixed, sorghum.power, sorghum.ConstPower, sorghum.exp)
##                    df      AIC
## sorghum.ols         3 18017.89
## sorghum.fixed       3 22316.61
## sorghum.power       4 17964.38
## sorghum.ConstPower  5 17966.38
## sorghum.exp         4 17921.13
#Exp
summary(sorghum.exp)
## Generalized least squares fit by REML
##   Model: logSorghum ~ rescale_ND 
##   Data: sorghum.near 
##        AIC      BIC    logLik
##   17921.13 17950.06 -8956.567
## 
## Variance function:
##  Structure: Exponential of variance covariate
##  Formula: ~rescale_ND 
##  Parameter estimates:
##       expon 
## -0.01997319 
## 
## Coefficients:
##                  Value   Std.Error   t-value p-value
## (Intercept) -0.5346044 0.012882963 -41.49701       0
## rescale_ND   0.0535188 0.001519261  35.22687       0
## 
##  Correlation: 
##            (Intr)
## rescale_ND -0.897
## 
## Standardized residuals:
##        Min         Q1        Med         Q3        Max 
## -4.9118332 -0.4762013  0.2568962  0.6495763  2.1113273 
## 
## Residual standard error: 0.6680175 
## Degrees of freedom: 10229 total; 10227 residual

#Soybean -

#Select Data and Rescale Distance
near<-read.csv(file="./Near_All/Soybean_Near_Centroid.csv", header=T)
near$rescale_ND <- near$NEAR_DIST/(1000*1000)
soybean.near<-near %>% select(mean_soybe,rescale_ND)
soybean.near<-soybean.near %>% filter(mean_soybe > 0, na.rm=TRUE)
soybean.near<-soybean.near %>% filter(rescale_ND > 0, na.rm=TRUE)
soybean.near$logSoybean<-log10(soybean.near$mean_soybe)

#Plot
ggplot(soybean.near, aes(x=rescale_ND, y=logSoybean)) +geom_point(alpha=0.2, color="#00A9B4") +geom_smooth(method = "lm",color="#DB3500") +theme_justin +xlab("Distance From GCO (1000's km)") +ylab ("Log10 Yield")
## `geom_smooth()` using formula 'y ~ x'
Soybean Near near Model,

Soybean Near near Model, 

###  Model - Soybean 
# regular OLS no variance structure
soybean.ols <- gls(logSoybean ~ rescale_ND, data = soybean.near)
# varFixed (variance changes linearly with X)
soybean.fixed <- update(soybean.ols, .~., weights = varFixed(~rescale_ND))
# varPower (variance changes as a power function with X)
soybean.power <- update(soybean.ols, . ~ ., weights = varPower(form = ~rescale_ND))
# varExp (variance changes as an exponential function of x)
soybean.exp <- update(soybean.ols, . ~., weights = varExp(form = ~ rescale_ND)) # error
# varConstPower (constant plus a power function of X (useful if X includes 0))
soybean.ConstPower <- update(soybean.ols, . ~., weights = varConstPower(form = ~ rescale_ND))
# compare all models by AIC
AIC(soybean.ols, soybean.fixed, soybean.power, soybean.exp, soybean.ConstPower)
##                    df      AIC
## soybean.ols         3 17503.72
## soybean.fixed       3 21913.61
## soybean.power       4 17499.79
## soybean.exp         4 17485.46
## soybean.ConstPower  5 17501.79
#varPower
summary(soybean.exp)
## Generalized least squares fit by REML
##   Model: logSoybean ~ rescale_ND 
##   Data: soybean.near 
##        AIC      BIC    logLik
##   17485.46 17514.52 -8738.731
## 
## Variance function:
##  Structure: Exponential of variance covariate
##  Formula: ~rescale_ND 
##  Parameter estimates:
##        expon 
## -0.006650274 
## 
## Coefficients:
##                  Value   Std.Error   t-value p-value
## (Intercept) -0.3391456 0.009434459 -35.94754       0
## rescale_ND   0.0164092 0.001066931  15.37981       0
## 
##  Correlation: 
##            (Intr)
## rescale_ND -0.821
## 
## Standardized residuals:
##        Min         Q1        Med         Q3        Max 
## -5.4419479 -0.4105329  0.3369870  0.6747373  1.6719120 
## 
## Residual standard error: 0.5796952 
## Degrees of freedom: 10548 total; 10546 residual

#Sunflower -

#Select Data and Rescale Distance
near<-read.csv(file="Near_All/Sunflower_Near_Centroid.csv", header=T)
near$rescale_ND <- near$NEAR_DIST/(1000*1000)
sunflower.near<-near %>% select(mean_sunfl,rescale_ND)
sunflower.near<-sunflower.near %>% filter(mean_sunfl > 0, na.rm=TRUE)
sunflower.near<-sunflower.near %>% filter(rescale_ND > 0, na.rm=TRUE)
sunflower.near$logSunflower<-log10(sunflower.near$mean_sunfl)

#Plot
ggplot(sunflower.near, aes(x=rescale_ND, y=logSunflower)) +geom_point(alpha=0.2, color="#FFB200") +geom_smooth(method = "lm",color="#DB3500")  +theme_justin +xlab("Distance From GCO (1000's km)") +ylab ("Log10 Yield")
## `geom_smooth()` using formula 'y ~ x'
Sunflower Near near Model,

Sunflower Near near Model, 

###  Model - Sunflower 
# regular OLS no variance structure
sunflower.ols <- gls(logSunflower ~ rescale_ND, data = sunflower.near)
# varFixed (variance changes linearly with X)
sunflower.fixed <- update(sunflower.ols, .~., weights = varFixed(~rescale_ND))
# varPower (variance changes as a power function with X)
sunflower.power <- update(sunflower.ols, . ~ ., weights = varPower(form = ~rescale_ND))
# varExp (variance changes as an exponential function of x)
sunflower.exp <- update(sunflower.ols, . ~., weights = varExp(form = ~ rescale_ND)) # error
# varConstPower (constant plus a power function of X (useful if X includes 0))
sunflower.ConstPower <- update(sunflower.ols, . ~., weights = varConstPower(form = ~ rescale_ND))
# compare all models by AIC
AIC(sunflower.ols, sunflower.fixed, sunflower.power,sunflower.ConstPower, sunflower.exp)
##                      df      AIC
## sunflower.ols         3 15907.51
## sunflower.fixed       3 19567.89
## sunflower.power       4 15829.34
## sunflower.ConstPower  5 15831.34
## sunflower.exp         4 15881.72
#varPower
summary(sunflower.ConstPower)
## Generalized least squares fit by REML
##   Model: logSunflower ~ rescale_ND 
##   Data: sunflower.near 
##        AIC      BIC   logLik
##   15831.34 15867.03 -7910.67
## 
## Variance function:
##  Structure: Constant plus power of variance covariate
##  Formula: ~rescale_ND 
##  Parameter estimates:
##        const        power 
## 5.417105e-06 8.663531e-02 
## 
## Coefficients:
##                   Value   Std.Error    t-value p-value
## (Intercept) -0.29574009 0.012994679 -22.758552       0
## rescale_ND  -0.00778374 0.001392125  -5.591268       0
## 
##  Correlation: 
##            (Intr)
## rescale_ND -0.894
## 
## Standardized residuals:
##        Min         Q1        Med         Q3        Max 
## -4.5510755 -0.4076590  0.2530632  0.7213569  1.6555183 
## 
## Residual standard error: 0.4766152 
## Degrees of freedom: 9302 total; 9300 residual

#Wheat -

#Select Data and Rescale Distance
near<-read.csv(file="Near_All/Wheat_Near_Centroid.csv", header=T)
near$rescale_ND <- near$NEAR_DIST/(1000*1000)
wheat.near<-near %>% select(mean_wheat,rescale_ND)
wheat.near<-wheat.near %>% filter(mean_wheat > 0, na.rm=TRUE)
wheat.near<-wheat.near %>% filter(rescale_ND> 0, na.rm=TRUE)
wheat.near$logWheat<-log10(wheat.near$mean_wheat)

#Plot
ggplot(wheat.near, aes(x=rescale_ND, y=logWheat)) +
  geom_point(alpha=0.2, color="#00A9B4") + geom_smooth(method = "lm",color="#DB3500")+theme_justin +xlab("Distance From GCO (1000's km)") +ylab ("Log10 Yield")
## `geom_smooth()` using formula 'y ~ x'
Wheat Near near Model,

Wheat Near near Model, 

###  Model - Wheat 
# regular OLS no variance structure
wheat.ols <- gls(logWheat ~ rescale_ND, data = wheat.near)
# varFixed (variance changes linearly with X)
wheat.fixed <- update(wheat.ols, .~., weights = varFixed(~rescale_ND))
# varPower (variance changes as a power function with X)
wheat.power <- update(wheat.ols, . ~ ., weights = varPower(form = ~rescale_ND))
# varExp (variance changes as an exponential function of x)
wheat.exp <- update(wheat.ols, . ~., weights = varExp(form = ~ rescale_ND)) # error
# varConstPower (constant plus a power function of X (useful if X includes 0))
wheat.ConstPower <- update(wheat.ols, . ~., weights = varConstPower(form = ~ rescale_ND))
# compare all models by AIC
AIC(wheat.ols, wheat.fixed, wheat.power, wheat.ConstPower, wheat.exp)
##                  df      AIC
## wheat.ols         3 22293.30
## wheat.fixed       3 27927.21
## wheat.power       4 22226.63
## wheat.ConstPower  5 22228.63
## wheat.exp         4 22226.87
#varPower
summary(wheat.power)
## Generalized least squares fit by REML
##   Model: logWheat ~ rescale_ND 
##   Data: wheat.near 
##        AIC      BIC    logLik
##   22226.63 22256.46 -11109.32
## 
## Variance function:
##  Structure: Power of variance covariate
##  Formula: ~rescale_ND 
##  Parameter estimates:
##       power 
## -0.07538861 
## 
## Coefficients:
##                   Value   Std.Error   t-value p-value
## (Intercept) -0.10244907 0.010269171 -9.976371       0
## rescale_ND   0.00872257 0.001297664  6.721749       0
## 
##  Correlation: 
##            (Intr)
## rescale_ND -0.87 
## 
## Standardized residuals:
##        Min         Q1        Med         Q3        Max 
## -4.7594172 -0.4459564  0.2490299  0.6773372  1.7349711 
## 
## Residual standard error: 0.6526273 
## Degrees of freedom: 12812 total; 12810 residual